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Abstract: The rare Higgs boson decay H —> Z 7 is forbidden at tree-level. In the Standard 
Model, it is loop-mediated through a W boson or a heavy quark. We analytically compute 
the QCD correction to the heavy quark loop, confirming earlier purely numerical results, 
that were obtained for on-shell renormalization. The small quark mass expansion of the 
decay matrix element contains only single-logarithmic contributions at each perturbative 
order, which is in contrast to the double logarithms observed in H —>• 77 . We investigate 
the numerical interplay of bottom and top quark contributions, and the dependence of the 
result on the renormalization scheme. 
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1 Introduction 

The discovery of the Higgs boson [1, 2] initiated a large research program aiming to de¬ 
termine the Higgs boson properties and to establish the details of the Higgs mechanism 
of electroweak symmetry breaking. Experimental measurements are up to now consistent 
with expectations from the Standard Model of particle physics. It is anticipated that the 
upcoming LHC data taking period at higher energy and luminosity will provide more pre¬ 
cise measurements and open up new observables that were previously inaccessible. In this 
context, the different decay modes of the Higgs boson play a very vital role. The decays to 
massive gauge bosons and to fermions are allowed at tree-level and provide direct measure¬ 
ments of the Higgs boson couplings. The rare decays H —> 77 and H —» Z 7 are forbidden 
at tree-level, they are mediated through loops containing massive particles [3-5]. As such, 
they are more sensitive to new physics effects from high energy scales than the tree-level 
dominated decay modes. 

The H —> 77 decay mode was among the most significant signatures in the Higgs boson 
discovery [1, 2], it has been measured in the meantime [ 6 , 7] with a relative precision of 
below twenty per cent. The branching ratio for H —> Z 7 , including the leptonic branching 
ratio of the Z boson, is considerably lower and only upper bounds could be established on 
it up to now [ 8 , 9]. Once established, this decay will provide access to a broader spectrum 
of observables than H —> 77 , since the decay of the Z boson to leptons will enable the 
study of spin-dependent particle correlations. It should be noted that the decay H -7 Z 7 
will be identified through a mass-cut on the pair of decay leptons, and that it should be 
considered to be a pseudo-observable [ 10 ]. 
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Higher order QCD corrections to H —> Z 7 from gluon exchange in the top quark loop 
were derived in [ 11 ] by performing a purely numerical evaluation of the relevant two-loop 
integrals in terms of five-dimensional Feynman parameter representations. The results 
derived in [11] use an on-shell renormalization for the top quark mass and the Yukawa 
coupling. Electroweak corrections to this decay are not known at present [10, 12]. It is the 
aim of this paper to rederive the QCD corrections to H —> Z 7 in an analytical form and 
to quantify uncertainties on them arising from scheme and scale dependence. 

Besides its phenomenological implications for H —> Z 7 , our calculation also provides 
an important subset of two-loop integrals relevant to the two-loop amplitudes gg —>• Hg 
and qg —>■ Hq with full top quark mass dependence. These amplitudes are known at 
present only at one loop [13], corresponding to the leading order in perturbation theory. 
For precision studies of the transverse momentum distribution of the Higgs boson and of 
Higgs-boson-plus-jet production, an effective held theory in the limit of infinite top quark 
mass is used commonly. In this approach, NLO QCD corrections were derived [14-17] 
and the calculation of the NNLO corrections is well-advanced [18-22]. The effective held 
theory description is however inappropriate at large transverse momenta, where the top 
quark loop is resolved by the recoiling jet; and it is precisely in this region that deviations 
from the Standard Model due to new heavy states could become visible. The calculation 
of NLO QCD corrections with exact top quark mass dependence is therefore recognized as 
high-priority aim [23, 24] , and the integrals derived here will be an important step towards 
it. 

This paper is structured as follows: in Section 2, we establish the notation and discuss 
the different contributions to the H -> Z 7 decay. Section 3 describes the calculation of the 
amplitude, including a detailed discussion of the relevant two-loop three-point integrals and 
of the renormalization. The numerical results are discussed in Section 4, and we conclude 
with Section 5. 


2 The H —> Z 7 decay in the Standard Model 


The Standard Model does not allow a tree-level coupling of the Higgs boson to a photon 
and a Z boson. The process 

H(q) -t Z(p 1 )'y(p 2 ) 

is mediated through a virtual particle loop, containing either a W boson or a massive 
quark [4, 5]. The Lorentz structure of its Feynman amplitude is constrained by gauge 
invariance to contain only a single scalar form factor: 

p/iv 

M = Aei^ipi, Ai) e 2 ,i/(p 2 , A 2 ) , ( 2 . 1 ) 

where the projector P flv is given by 


i^=p£rf-(pi-P2)«r. 

The decay width H —> Z 7 is obtained as 


r = 


Gi a m 


2 

w 


4 m? H (rnjj — rn^) 




( 2 . 2 ) 


(2.3) 
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with Fermi’s coupling constant Gp, the fine-structure constant a and the Higgs and Z boson 
masses mp and m z , respectively. Depending on the particle coupled to the external Higgs 
boson, the amplitude can be further decomposed into contributions from the W boson and 
the fermions q: 

A = C\vA\y + c q A q ■ ( 2 - 4 ) 

q 

The coupling factors are 

2 Qq {Jq ~ 2 Qq sin2 ®w) 


C\y = cos 9 W , 


C q = N r 


(2.5) 


cos 9 W 

where 9 W is the weak mixing angle, N c the number of colors, Q q the charge of the fermion 
and Iq the third component of its weak isospin. Due to the mass hierarchy of the particles 
involved, we will only consider the dominant pieces coming from the W boson, the top 
quark and the bottom quark. 

The Born-level contribution to the amplitude arises at one loop, it is written as: 


n) 1 ) — cw aL) + ctA+ ^ + CbA 


(i) 


t(i) 


( 2 . 6 ) 


Higher-order perturbative corrections are obtained by a loop expansion of the amplitude A. 
Next-to-leading order QCD corrections affect only A t and Ab, they correspond to two-loop 
graphs with an internal mass: 

Aq(m H ,m Zl m qi a s , n) = AA\m H , m z , m q ) + -— M - A ( 2) (m H ,m z ,m q ,fj,). (2.7) 


( 2 ) 

For the dominant top quark contribution, the two-loop correction A\ ' was computed nu¬ 
merically (based on the Feynman parameter representation of the amplitude) using an 
on-shell renormalization for the quark mass and the Yukawa coupling in [11], 


3 Calculation of the amplitude 

To compute the one-loop amplitudes A^J and A ^ as well as the two-loop QCD contri- 
bution A q ' to the quark-mediated amplitude for H —> we project all relevant Feyn¬ 
man diagrams generated by Qgraf [25] onto the tensor structure (2.1) using Form [26]. 
The resulting Feynman integrals are reduced to a set of master integrals with the help of 
integration-by-parts (IBP) identities [27], which are solved using the Laporta algorithm [28] 
implemented in the Reduze2 code [29, 30]. After the reduction, the amplitude can be ex¬ 
pressed in terms of a certain number of master integrals depending on the loop order. 

Each of these master integrals has a specific mass dimension, which can be scaled out 
by multiplying with the appropriate power of the mass mi running in the loop, such that 
the resulting dimensionless integrals are only functions of the mass ratios m 2 H /mf and 
m? z /rnf . We parametrize this dependence by introducing Landau-type variables 


2 2(1 hiY 

m H = -rm 


hi 


2 2(1 z iY 

= ~mi - 


(i = W,q) 


(3.1) 


For the sake of readability, we drop the subscript q of these variables in the following 
whenever we deal with quark-mediated amplitudes, i.e. h q = h, z q = 2 . 
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Figure 1. Example diagrams (a) for the computation of the one-loop amplitudes A 1 ^) and Aq 1 ^ 
and (6) for the computation of the two-loop amplitude A\ . 


3.1 The one-loop amplitude 

The Born-level one-loop amplitudes for the contributions of IT-bosons A^} and heavy 
quarks Aq to H —>• Z'y were derived in [4, 5]. Example diagrams are depicted in Fig. 1(a). 
In terms of the Landau variables introduced above, these results read 


A 


(i) 4 % S e m^y 

W h W Zyy 


{(^ + l)(^ + l)-4 h w (z w + lf} | (hw ~ z w ){h w z w - f) 

~ ^ + hw- 1 — — log ( hw ^ + ( z w ~ lo s( 2 w)| 


q 


+ {(zw + Z\v( z W + 4))(/ipp + 1) — 2 hw(zw — l) 2 } 
{log 2 (h w ) - log 2 ( 2 p/)}] + 0 (e), 


i S e m 3 y q v 


1 


1 


8Vl+ h^{ Z + ~z 

(h + l)(z- l) 2 ( 1\ 

' z{h — 1) lo gW+^-^j lo ^) 

- 2 {li+^-^--+ 4 } (log 2 (/j.) - log 2 (z)} 


+ 0 (e) 


(3.2) 


(3.3) 


with the Standard Model Higgs vacuum expectation value v and the Yukawa coupling 


m n 


Vq = — 
v 

associated with the quark q. The normalization factor 

?T (1 + e) ( 47T/Xo V 


S f = 


167r 2 \ m 2 J 


(3.4) 


(3.5) 


arises from the integration measure f d D k/(2n) D of the master integrals in D = 4 — 2e 
dimensions, where fi$ is the mass scale of dimensional regularization. 


3.2 Differential equations and integral basis 

( 2 ) 

The two-loop amplitude for the quark contribution A q ; has been computed purely numer¬ 
ically in terms of a five-dimensional Feynman parameter integral in [11], We derive an 
analytical expression for this amplitude, through a reduction of all two-loop integrals to a 
set of master integrals. 
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To compute the two-loop master integrals, we use the method of differential equa¬ 
tions [31-34], In this method, differential equations in internal masses and external in¬ 
variants are derived for each integral by performing the differentiation on the integrand, 
which is then related to the original master integral by the IBP identities. With this, 
we obtain inhomogeneous differential equations in either Landau variable, plus a trivial 
homogeneous equation in m q for each integral. The differential equations are solved in a 
bottom-up approach, i.e. starting from the master integrals with the lowest number of 
different propagators (‘topology’) because they will show up in differential equations of 
higher topologies. 

The coefficients of the individual master integrals in the homogeneous and inhomoge¬ 
neous terms of the differential equations are rational functions of h and 2 . Upon partial 
fractioning, only a limited number of polynomials in h and z appear. These form the 
so-called alphabet associated with this set of master integrals: 

12 } (3.6) 

with 

h = z, 

I 2 = z T 1, 
h = z-1 , 

14 = h, 

15 = h + 1, 

1% — h 15 
£7 = h — z , 

Is = hz — 1, 

£9 = h 2 — hz - 

£10 = h 2 z — hz 
In = z 2 — hz - 
I 12 = z 2 h — hz 

With an appropriate choice of basis integrals [34] , the full system of differential equa¬ 
tions for all 28 master integrals (written as 28-component vector M) takes the form of a 
total differential, 

12 

d M(h,z) = ey^ R k d\og(l k ) M(h, z ), (3.7) 

k =1 

where the matrices R k contain only rational numbers. In this case, two important features 
can be exploited. First, the differential equations can be integrated order by order in e 
in terms of generalized harmonic polylogarithnrs (GHPLs) [35-38]. They are defined as 


■ h + 1, 

— h + z, 

— z + 1, 

— z + h. 
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iterated integrals according to 


f x 1 

G(wi,...,w n -x)= / d t- - G(w 2 , ■ ■ ■ ,w n ;t) , 

Jo t ~ w i 

log n x 


G 


^ 0 n ; xj 


n\ 


(3.8) 


with indices Wi and argument x. Second, the results will be expressed as a linear com¬ 
bination of GHPLs of homogeneous weight. We arrive at a total differential of the form 
(3.7) starting from the Laporta basis / depicted in Fig. 2 and subsequently applying the 
algorithm described in [39]. The resulting canonical basis M in terms of the Laporta basis 
/ is given in Appendix A. 

Since the alphabet (3.6) is not linear in the Landau variables, we further decompose 
it to enable the integration in either h or z, which yields a solution up to an integration 
constant that only depends on the other variable. This boundary value is then determined 
by imposing regularity in special kinematic points, where the integrals are known to be 
regular from physical arguments. In our case, these points are given by 


h = 

1 

■fA 


= 0, 


z = 

1 

•H- 

m 2 z 

= 0, 


h = 

z 

■H- 


= m 2 z , 


h = 

1 

‘V 

■H- 


= rn 2 z , 

(3.9) 


i.e. they correspond to the limit where the masses of the external particles either vanish 
or coincide. Choosing one of these points such that the rational prefactors in Eqs. (A.l) 
are equal to zero considerably reduces the complexity of the integration constants. 

By taking limits in these kinematic points, we are left with GHPLs that contain the 
same variable x G {h, z} both in the argument and in the indices. In order to simplify the 
result and to obtain a unique representation, we use an inhouse MATHEMATICA implemen¬ 
tation [40, 41] relying on the symbol and coproduct formalism [42] to transform GHPLs of 
the type 

G (u>i(x),... ,w n (x)-x) -> G (ai,..., a n \ x) , (3.10) 

where the transformed indices ai are nothing but complex numbers. In doing so, we end 
up with GHPLs up to weight four, which are given by 


G (cii,..., a ni L) 
G (bi,... ,b n ; z) 


with ai G {0, ±1, K ±, Lf } , 

Z Jz 

with bi G {0, ±1, c, c} , 


(3.11) 
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Figure 2. Two-loop master integrals for the calculation of Aq*K Dashed lines are massless, whereas 
internal solid lines denote propagators with mass m q . Double and solid external lines correspond 
to virtualities rn 2 H and respectively. Dotted propagators are taken to be squared. 
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where 


1 / 

c = j (l + l ^V ’ 

j = _ z _ 

z 1-z + z 2 ’ 

Kf = \{^ + z± y/-3 + 2 z + z 2 ^j , 

Lf = ^ (l + ^ ± \J\ + 2 z - 3 z 2 ') (3.12) 

for the underlying set of master integrals. Further transformations of the type 

G(w 1 (x),.. .,w n {x)\y) -» G{ci(y), ... ,c n (y);x) (3.13) 

become necessary when the integration is performed in a different variable compared to 
the integration of a master integral of a subtopology which enters the differential equation 
under consideration. 

It remains to comment on one issue: Transformations of the types (3.10) and (3.13) 
are performed for all master integrals except for M 25 —M 2 S, where GHPLs of the form 

G(wi(x),...,w n (x)-,x) (3-14) 

with Wi = {K±,L±} occur. The transformations were not necessary in this case because 
the results of these integrals do not enter the differential equation of any master inte¬ 
gral of higher topology. However, transformations of these types could become desirable 
when four-point functions are computed, where these four coupled master integrals appear 
as subtopologies. This could be attempted by a multivariate extension of the general¬ 
ized weights approach described in [43], i.e. rather working on non-linear indices of the 
form (3.6) than on linear ones as in Eq. (3.11). 

We would like to state that the results of the master integrals were checked in several 
ways. We performed transformations of the type (3.13) and verified that the solution fulfills 
the differential equation in the other variable. This check works only up to a constant, 
which is why we compared each master integral numerically against SecDec [44] and 
found agreement to high precision. The analytic expressions of the master integrals are 
rather lengthy and will not be reproduced here. They are available in Mathematic A and 
Form format together with the arXiv submission of this paper. 

3.3 Calculation of the two-loop amplitude 

For the two-loop amplitude, six generic Feynman diagrams and their permutations have to 
be evaluated. They emerge from the one-loop diagram in Fig. 1(a) by attaching one gluon 
propagator to the fermion lines in every possible way, which is depicted in Fig. 1(6). After 
the manipulations described in the beginning of Section 3 and after inserting the analytic 
results of the master integrals, we are left with the unrenormalized two-loop amplitude. 

For its renormalization, we consider three different prescriptions: 

(a) quark mass and Yukawa coupling in the on-shell (OS) scheme. 
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(b) quark mass in the OS scheme, Yukawa coupling in the MS scheme. 


(c) quark mass and Yukawa coupling in the MS scheme. 


Since the two-loop amplitude is the leading order in a s , no renormalization of the gauge 
coupling is required. 

All three prescriptions yield the same pole parts of the renormalization counter terms 
and produce finite expressions for the renormalized amplitude. They are related by finite 
scheme transformations, which is why we choose to compute the renormalized amplitude in 
scheme (a) and use it to derive the results in schemes (b) and (c). In the pure OS scheme, 
the quantity 

— SmosCW + ZosAW (3.15) 

m q H H 

has to be added to the unrenormalized two-loop amplitude in order to remove its diver¬ 
gences, where Zos and <5mos are the Yukawa and mass renormalization constants, respec¬ 
tively [45]: 


Zos 

6m 0 s 


Osip) 

7 r 


16 i 7r 2 S e 


Cf 

4 


3 - 2e 

e (1 — 2e) 


m q Z os • 


(3.16) 

(3.17) 


Note that Eq. (3.15) requires the calculation of the one-loop amplitude A ^ and the mass 
counterterm C q up to O(e). Cq 1 ' 1 can be computed from the sum of the three diagrams 
that are obtained by putting a mass insertion into one of the fermion lines in Fig. 1(a). 

Next, we express the OS quantities M q and Y q in terms of MS quantities fn q and y q at 
a particular matching scale y rn using the standard relations (e.g. [46-48]): 


M q = m q (n) (1 + A) , 
Yq = Vq{v) (! + A ) , 


A = 


a s (ii) 


7T 


C F 


3 y 

1 + 7 lo g =2~ 




(3.18) 


We perform this matching at the scale of the running MS quark mass m q . Starting from 
the OS result A^’ a (m#, mz, M q ), these scheme transformations induce finite shifts in the 
amplitudes of the prescriptions (b) and (c): 


A^’ b \m H ,m z ,m q ,fi) = A ( q 2 ’ a) (m H , rn z , m q (fi)) + A • A^\m H , m z , m q (n )), 
A^ c \mH,m z ,m q ,n) = A^’ b \m H ,m z ,m q (y)) + A - ° Aq 


(3.19) 

M q =m q (n) 


In practice, the coefficient A emerges by making the following replacements in Eq. (3.3), 
where the Landau variables h and z are defined according to Eq. (3.1) with rn q = m q (fj,): 

h = h — 2 A h =-, 

h + 1 
z — 1 

z = z — 2Az -. (3.20) 

z + 1 v ; 
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The amplitudes in the schemes (a), (b) and (c) have a common polynomial structure, 
which contains only a limited number of combinations of the denominators (3.6): 



16-7T 2 S' 2 rn q y q v Cp • 


Cl , C 2 1 C3 | C4 i C5 | Cq | C7 | C8 

h + JJ ~ 5 + TJ 6 + hM + hi 9 i w + hh + hh + k 

C9 C10 Cn C12 C13 Cl4 Cl5 Cl6 

Ms M11 hh hh h h ho h2_ 


+ 0(e). 


(3.21) 


The coefficients c, are linear combinations of GHPLs multiplied by some power of h or 
z with positive exponent. The complete analytic expression of the two-loop amplitude 
exceeds the scope of this paper and is attached to the arXiv submission of this article. 

The analytic result of the full two-loop amplitude enables us to derive its limit for 
small quark masses and gain information about its logarithmic structure. We perform this 
expansion by removing trailing zeros in the indices of the GHPLs using the shuffle relation 
so that the logarithmic singularities become explicit. Subsequently, we apply the scaling 
relation to the remaining GHPLs of the form (3.11), which turn into 

1 r 1 z 1 J z 1 Kt 

G a„4) With 

G (hi,- b n ] 1) with bi € {0, ±-, 4 -} , (3.22) 

where a n ^ 0 and b n ^ 0. This procedure shifts the dependence on the quark mass from the 
argument to the indices of the GHPLs and allows expanding the integrand of their integral 
representation without the need to take care of the integration boundaries. Consequently, 
we solve the definition of the Landau variables in Eq. (3.1) for h and z, replace them in 
the integral representation and expand the result in m q . 

A subtlety occurs when GHPLs of the form (3.14) are analyzed for small quark masses. 
The particular case w\ = L~ gives rise to further singularities due to 


lim = z + 0(z 2 ), 

m q —hO 


(3.23) 


which is why the corresponding GHPLs are isolated with the help of the shuffle relation. 
Finally, we expand the remaining GHPLs with Wi = {/i^,L^}(i ^ 1) in small quark 
masses, apply transformations of the type (3.10) and treat the resulting GHPLs in the same 
way as the ones above. In doing so, we obtain the following expression for the amplitude 
from Eq. (2.7) in the limit of small quark masses, renormalized in the OS scheme: 


( M 2 \ ( m? \ 

—f- log I —f ) (3.24) 

) \ m H / 


1 + c F 2M log h$ 

vr \m 2 z 


Wlr 


4-8 l0g l^J + 2 l0 H 1_ ^ ,+ 


1717 


Ll2 


~ C 2 


2 log 
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To check the expansions, we validated numerically (using the GiNaC [38, 51] implemen¬ 
tation) that the individual GHPLs converge towards their expansions in the limit of small 
quark masses. In addition, we rederived the corresponding limit of the H —> 77 amplitude 
by starting from the H —> Z 7 amplitude and setting the Z boson mass to zero. With 
the manipulations described above, we agree with the previously available ratio of the 
two-to-one-loop amplitude for the process H —> 77 in small quark masses [13, 49]. 

It is important to note that Eq. (3.24) contains only single-logarithmic terms, which 
is in contrast to the H —> 77 case. The non-trivial cancellation of the double-logarithmic 
terms, which originate from the Sudakov region, takes place both in the one-loop and in 
the two-loop amplitudes of Eqs. (3.3) and (3.21), leaving only single-logarithmic terms at 
each order. A double-logarithmic Sudakov resummation, as performed for H —> 77 in [49] 
is therefore not needed for H —> Z 7 . We observe from Eq. (3.24) that the introduction of 
a running Yukawa coupling, scheme (b), resums the single-logarithmic contribution that is 
independent of m 2 z /m 2 H . 

4 Numerical results 

The calculation of the master integrals and the amplitude outlined in Section 3 is performed 
for the values 

0<h<l, 0 < z < 1. (4.1) 

This corresponds to the Euclidean region, where m 2 H and m 2 z in Eq. (3.1) are negative 
and the master integrals are real. In order to get a physical expression, the results have 
to be analytically continued to the physical Minkowski region, where we distinguish three 
kinematic regions: 

• Region I: m 2 z < m 2 H < 4 m 2 , 

• Region II: m 2 z < 4 m 2 < m 2 H , 

• Region III: 4 m 2 < rn z < m 2 H . 

In Region I, the virtualities of the external massive particles are below the threshold induced 
by the particle running in the loop and the amplitude is real [50]. The two solutions of 
each variable in Eq. (3.1) become imaginary in this region with 

h = e ± 2 i$ H ^ 

z = e ±2i * z . 


(4.2) 

(4.3) 


4>h and 4>z are phase factors given by 


4>i = arctan 



(4.4) 


i.e. the variables lie on the unit circle in the complex plane and the second solution is the 
complex conjugated of the first one. In general, care has to be taken when choosing one of 
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them such that it is in agreement with the common +i0 prescription for the Mandelstam 
variables m 2 H and m 2 z . Due to the real amplitude, however, there is no ambiguity in this 
case and the imaginary parts of h and z can be chosen freely. 

In Region III, w? H and m 2 z are beyond the threshold and the amplitude picks up an 
imaginary part. The two solutions of each variable in Eq. (3.1) are real, with one fulfilling 

— 1 < h < 0 , -l<z<0, (4.5) 


and the other one being its inverse. In contrast to Region I, one has to be careful when 
assigning a small imaginary part to h and z in order to fix branch cut ambiguities related 
to the Mandelstam variables. In this case, a positive imaginary part leads to the correct 
result if \h\ < 1 or \z\ < 1 and a negative imaginary part has to be applied when \h\ > 1 or 
\z\ > 1. 

From the input parameters specified below, it is obvious that the top quark amplitude 
is calculated in Region I, while the bottom quark amplitude is computed in Region III. 

Region II is not needed for the physical values of the masses. 

( 2 ) 

Since the analytical expression for A q is given in terms of GHPLs, it can be evaluated 
using GiNaC [38, 51]. For masses and couplings, we use the input values 


af\m H ) = 0.1130114, 
rriH = 125.7 GeV , 
M t = 173.21 GeV, 
hh,(m#) = 2.7832 GeV, 
Qt = 2/3, 

4 3 = - 1 / 2 , 


a = 1/128, 
m z = 91.1876 GeV, 
m((mjj) = 167.21 GeV , 
sin 2 6 W = 0.23126 , 

Qb = ~ 1/3, 

C F = 4/3, 


G f = 1.1663787 • 10" 5 GeV" 2 , 
mw = 80.385 GeV , 

M b = 4.7652 GeV , 
cos 2 6 W = 0.76874, 

It= 1 / 2 , 

N c = 3. (4.6) 


They were obtained by evolving the values of the Particle Data Group Collaboration [52] 
to the scale /r = with the two-loop renormalization group equations [53]. To resum 
potentially large single logarithms in m q /m.H to all orders in the perturbative expansion, 
we use the two-loop renormalization group equations [53] to evolve the MS quark mass 
(and accordingly the Yukawa coupling) from the matching scale to /r = rrijj. This leads 
to the following next-to-leading-order decay width in the renormalization schemes (a), 
(b) and (c): 


p(2,a) 


7.075:1! + 0.42800 




7r 


keV 


M = H 7.09072keV. 


p(2,6) 


7.09409 + 


a s {m H ) 


TV 


-0.53266 - 0.76661 log — 


m 2 (m/f) 


+ 0.01229 log — 


= 7.09403 keV. 


(4.7) 

_' 2 ) 
m 2 {m H ) J _ 

(4.8) 


keV 


p(2,c) 


7.05934 + 


a s {m H ) 


7T 


0.64587 + 0.10597 log 


m 


H 


m 2 {mH) 


+ 0.01453 log 


m 


H 


= 7.08438 keV. 


m 2 {m H ) 

(4.9) 


keV 
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The breakup of the terms in the first lines of Eqs. (4.8) and (4.9) is to illustrate the relative 
numerical importance of the individual contributions. 

An estimation of the uncertainty on the prediction from missing higher orders is pro¬ 
vided by varying mu /2 < fi < 2 mu in Fig. 3. For every data point /i = fj,Q, the MS quark 
mass rn q (fi o) and the strong coupling constant a s (/j, o) are evolved to the scale no using 
the two-loop renormalization group equations [53]. The leading-order decay width T^ 1 ) 
and the next-to-leading-order decay width r^) can be separated into contributions from 
the W boson, top quark and bottom quark amplitudes as well as their interferences. The 
corresponding values are shown in Table 1, from which it becomes clear that the bottom 
quark amplitude has to be taken into account, since its interference with the W amplitude 
is of the same order of magnitude as the self-interference of the top quark amplitude at 
one loop. Moreover, the combination T^/ fe exceeds the top quark self-interference T) f at 
two loops. 

Our on-shell results are in agreement with the numerical findings of [11], Furthermore, 
we performed a detailed numerical comparison with [54] and found agreement to high 
precision. 


Table 1. Various contributions to the numerical result of the leading-order decay width Ob and the 
next-to-leading-order decay width O 2 ) in the renormalization schemes (a), (b) and (c), evaluated for 
/r = to j-i . In case of I']b, the subscripts i and j indicate the interference of the one-loop amplitudes 

Ayy , A.] 1 ' and A^\ whereas F <2> describes the interference of the one-loop amplitude with the 
( 2 ) 

two-loop amplitude A) '. All values are given in keV. 


Partial width 

(a) 

(b) 

(c) 

p(l) 

1 ww 

7.86845996 

7.86845996 

7.86845996 

p(t) 

1 Wt 

-0.83636436 

-0.80736905 

-0.84015333 

p(l) 

1 Wb 

0.02216139 

0.01294390 

0.00908488 

p(l) 

1 tt 

0.02222498 

0.02071068 

0.02242680 

p(t) 

1 tb 

-0.001177803 

-0.00066408 

-0.00048502 

p(t) 

1 bb 

0.00002103 

0.00000717 

0.00000325 

r<4) 

7.07532519 

7.09408860 

7.05933655 

p(2) 

1 Wt 

0.02213199 

-0.00078617 

0.02467587 

p( 2 ) 

1 Wb 

-0.00588750 

0.00073044 

0.00176120 

p(2) 

1 tt 

-0.00117624 

0.00004033 

-0.00131738 

p(2) 

1 tb 

0.00031290 

-0.00003747 

-0.00009403 

p(2) 

1 bt 

0.00003117 

-0.00000065 

0.00001425 

p( 2 ) 

1 bb 

-0.00001592 

-0.00000081 

0.00000078 

r( 2 ) 

7.09072159 

7.09403427 

7.08437723 
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(i in GeV 

Figure 3. Scale variation of the next-to-leading-order decay width C 2 ) in the renormalization 
schemes (a), (b) and (c) for 30 GeV < /j, < 300 GeV. 


We observe that the NLO results for the decay width are consistent between the three 
schemes. The relative size of the NLO correction is 2 % 0 in scheme (a), below ICY ' 5 in 
scheme (b) and 3%o in scheme (c). The very small corrections in scheme (b) are however 
in large part due to numerical cancellations between a priori unrelated terms. The spread 
between the different schemes is 1.3% 0 at fi = mu , and variations of the renormalization 
scale change the predictions in either given scheme by at most 0.4% o . 

5 Conclusions 

In this paper, we have revisited the QCD corrections to the rare loop-induced Higgs boson 
decay H —> Z'y. The relevant two-loop three-point integrals with two different external 
masses and one internal mass were derived analytically, using a reduction to master in¬ 
tegrals, which were then computed using differential equations. These integrals are also 
an important ingredient to the two-loop amplitudes for Higgs-plus-jet production in gluon 
fusion with full dependence on the internal quark masses. 

By expanding the one-loop and two-loop matrix elements in the on-shell scheme for 
H —> Z 7 in the limit of small quark masses, we noted the absence of double-logarithmic 
contributions, which is in contrast to H —> 77 [13, 49]; single-logarithmic terms are re¬ 
summed by the introduction of a running Yukawa coupling. 
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We investigated the dependence of the corrections on the renormalization scheme used 
for the quark mass and Yukawa coupling. We observe that the results for the decay rate in 
on-shell and MS schemes, as well as in a hybrid scheme with on-shell mass and MS Yukawa 
coupling, are well consistent with each other, and that corrections are in the sub-per-cent 
range in all three schemes (being smallest in the hybrid scheme). We confirm the previously 
available numerical on-shell result [11] and agree with an independent caclulation [54], The 
residual QCD uncertainty on the H —> Z 7 decay rate is around 1.7%o from the combination 
of scale variation and spread between the different renormalization schemes. 
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A Canonical Master integrals 


In this appendix, we provide the relations between the two-loop canonical master integrals 
that appear in Eq. (3.7) and the Laporta master integrals of Fig. 2 in terms of the Landau 
variables defined in Eq. (3.1). The canonical integrals are normalized such that their 
Laurent expansion starts at order e°. In addition, we extract the mass dimension of the 
integrals and obtain 


Mi 


M 2 


m 3 


m 4 


m 5 


m 6 


m 7 


M s 


e 2 h, 


2 2 

e rn q 


2 2 

e rn q 


2 2 
e m q 


2 2 

e rn q 


2 2 

e rn q 


2 2 
e m q 


2 4 
e m q 


(z + l)(z-l) T 

-A, 

z 

— — — [( z + 1 ) -£3 + -£4] > 
(* - l) 2 j 

(h + 1 )(h — 1) 

- h - 

—^— i(h + 1 ) h + h\ > 

(h- l) 2 

- h - 17 ’ 

(z + l) 2 (z-l) 2 r 

- 2 - 7 8 > 
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Mg = — e 3 m 
Mio = — e 3 m 
Mn = e 2 w-J ■ 
M 12 = e 2 m 4 ■ 
M13 = —e 3 ?n 

M m = — e 3 m 
A/i 5 = —e 2 m 

Afig = —e 3 m 
Mn = —e 3 m 
M\s = —e 2 fn 

Mn = e 3 m 4 
M 2 o = -e 3 m 
M 21 = — e 4 m 
M 22 = -e 3 m 
M 23 = -e 4 m 
M 21 = —e 3 m 
M 25 = e 4 m 2 ■ 
M 26 = e 3 m 4 ■ 
M 27 = e 3 m 4 
M 28 = e 2 m 4 ■ 


2 ~~ z)(hz — 1 ) 

q hz 19 ’ 

2 (/»-«)(/»«- 1 ) r 
» hz /l0 ' 

(fr 2 -l)(z 2 -l) 
hz 

(h + l) 2 (h-l) 2 
h 2 

2 0 - z)(hz - 1 ) 

9 /iz 

2 (h - z)(hz - 1 ) 
q hz 

2 z ( h 2 + 1 ) — h (z + 1 ) 

9 2z(h 2 + l)~h(z + l) 2 


2 (h-z)(hz- 1) r 

9 - hz - 16 ’ 

2 (h — z)(hz — 1 ) 

- hz /l7 ’ 

2 h (z 2 + 1 ) — z (h , + 1 ) 

9 2 h(z 2 + 1 ) - z {h + l ) 2 


Al 5 

A 2 , 

^13 : 
Jl4 . 


3 (h — l ) 2 (h — z)(hz — 1 

2 s ^^ 

(z 2 -1)(/7 2 + 1-/i(z + 1)) 


+m; 


/iz 


3 (z — 1)“ (h — z)(hz — 1 / 

- 74 — e---- 

2 z faz 

(/i 2 -l)(z 2 + l-z(/i + l)) 


+mt 


/iz 


(/i — z)(hz — l)(z 2 — 1 ) 
hz 2 

4 (h — z)(hz — 1 )(h 2 — 1 ) 
9 Wz 

2 {h - z)(hz - 1 ) 

9 hz 


119 , 
-^20 


-^21 


(h — z)(hz — 1 )(h 2 — 1 ) 


9 h 2 z 

2 (h - z)(hz - 1 ) 


^23 


9 /lZ 

4 (h - z)(hz - l)(z 2 - 1 ) 
9 ^2 

(h — z)(hz — 1 ) 

T -*25 5 

hz 

(h- z)(hz- l)(z 2 - 1 ) 
hz 2 

(h — z)(/iz — l)(/i 2 — 1 ) 


1-22 


I2A , 


h 2 z 


1 26 , 

1 27 , 


hz — 1 
hz 


—4 (hz — 1) In + 2 e (/i — z) 


z - 1 


-^26 — 


/»- 1 


^27 


- (2 il3 + -^ 14 ) 

-^15 j. 


- (2 /i 6 + -^17) 

-^18 j 
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(A.l) 


+m 


2 (h - z) 2 {hz - 1) 
hz 


hs 


I± is the two-loop tadpole integral with both propagators taken to be squared. It is given 
by 


d D k f d u l 


h = 


D, 


1 


so that 


( 2 ^)-° J (2ir) D ( k 2 — m 2 ) 2 ( l 2 — to 2 ) 2 e 2 


M 1 = S 2 


(A.2) 


(A.3) 


where S 2 is the common normalization factor of all two-loop master integrals defined in 
Eq. (3.5). 
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